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ABSTRACT 

Recent observations have shown that star-forming galaxies like our own Milky Way evolve kinemat- 
ically into ordered thin disks over the last ~ 8 billion years since z = 1.2, undergoing a process of 
“disk settling.” For the first time, we study the kinematic evolution of a suite of four state of the art 
“zoom in” hydrodynamic simulations of galaxy formation and evolution in a fully cosmological con- 
text and compare with these observations. Until now, robust measurements of the internal kinematics 
of simulated galaxies were lacking as the simulations suffered from low resolution, overproduction of 
stars, and overly massive bulges. The current generation of simulations has made great progress in 
overcoming these difficulties and is ready for a kinematic analysis. We show that simulated galaxies 
follow the same kinematic trends as real galaxies: they progressively decrease in disordered motions 
(cr ff ) and increase in ordered rotation (Wot) with time. The slopes of the relations between both a g 
and V ro t with redshift are consistent between the simulations and the observations. In addition, the 
morphologies of the simulated galaxies become less disturbed with time, also consistent with observa- 
tions. This match between the simulated and observed trends is a significant success for the current 
generation of simulations, and a first step in determining the physical processes behind disk settling. 
Subject headings: galaxies - formation, galaxies - evolution, galaxies - kinematics and dynamics, 
galaxies - fundamental properties 


1. INTRODUCTION 

Over the last ~ 8 billion years since a redshift of one, 
the population of star-forming galaxies of ^Milky Way 
mass has settled kinematically into flat, rotationally- 
supported, disk galaxies (Kassin et al. 2007, 2012). In 
the past, these galaxies had more disordered motions (as 
measured via an integrated gas velocity dispersion, a g ) 
and less ordered rotation (V rot ) than they do today (Flo- 
res et al. 2006; Kassin et al. 2007; Vergani et al. 2012; 
Kassin et al. 2012). Over 0.1 < z < 1.2, the median a g 
of star-forming galaxies of ^Milky Way mass has pro- 
gressively decreased while V rot has increased (Kassin et 
al. 2012). Observations at even higher redshifts (z ~ 
1.5 — 3), albeit where galaxy samples are much less rep- 
resentative, also find star-forming galaxies with large 
amounts of disordered motions, as measured through 
<j g (e.g., Forster-Schreiber et al. 2009; Law et al. 2009; 
Wright et al. 2009; Lemoine-Busserolle & Lamareille 
2010; Lemoine-Busserolle et al. 2010; Gnerucci et al. 
2011 ). 

These findings are an important benchmark for any 
theory or simulation of disk galaxy formation. Con- 
straints on the internal kinematics of galaxies are sig- 
nificantly more stringent than constraints on e.g., lu- 
minosity, stellar mass, or star formation rate. Inter- 
nal kinematics directly dictate how the gas in galaxies 
is arranged, which is heavily influenced by the physical 
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processes involved in galaxy evolution. Furthermore, if 
theory is able to successfully reproduce observations of 
kinematic evolution, then it becomes a useful tool to help 
interpret the observations. In particular, we would like 
to know the physical processes behind disk settling. 

Until now, robust measurements of the internal kine- 
matics of simulated galaxies were lacking as the simu- 
lations suffered from low resolution, overproduction of 
stars, and overly massive bulges. These problems were 
ameliorated by increased resolution (e.g., Mayer et al. 
2001), better treatment of feedback (Stinson et al. 2006; 
Hopkins et al. 2012; Wise et al. 2012), and modeling of 
the ultraviolet background (Quinn, Katz, & Efstathiou 
1996). The current generation of galaxy simulations is 
able to match the rotation velocities, sizes, stellar masses, 
and baryon fractions of local disk galaxies (e.g., Gover- 
nato et al. 2007; Brooks et al. 2011; Brook et al. 2012; 
Munshi et al. 2013; Christensen et al. 2013; Vogelsberger 
et al. 2013), as well as some properties of higher redshift 
galaxies (e.g., Hirschmann et al. 2013; Zemp et al. 2013; 
Hopkins et al. 2013). They are the first generation of 
simulations which are able to make quantitative predic- 
tions for the internal kinematics of galaxies. 

However, most simulations have concentrated on repro- 
ducing only the rotation velocities of local disk galaxies 
(e.g., Guedes et al. 2011; Stinson et al. 2013; McCarthy 
et al. 2012; Brook et al. 2012). A few have measured a g , 
but only for z = 0 (e.g., Roskar et al. 2013; Agertz et 
al. 2013) or z = 2 (e.g., Ceverino, Dekel, & Bournaud 
2010; Angles- Alcazar et al. 2013). Croft et al. (2009) 
measured V rot and a g at z = 0 and z = 1 for galaxies in 
their simulations, and made comparisons to the observed 
V rot and So . 5 = y/0.5V r „ t + <r 2 Tully-Fisher relations from 
Kassin et al. (2007). No studies have yet compared sim- 
ulated galaxies to the progressive disk settling found at 
0.1 < z < 1.2 by Kassin et al. (2012). 
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In this paper, for the first time, we determine whether 
fully cosmological simulations of star-forming galaxies 
undergo the kinematic disk settling found in observations 
of real galaxies. 

2. HYDRODYNAMIC SIMULATIONS 

The kinematics of star-forming galaxies in the real uni- 
verse are measured from nebular emission lines which 
trace gas heated in star-forming regions. In order to 
compare hydrodynamic simulations of galaxies most di- 
rectly with observations, the simulations need to make 
predictions for the gas in galaxies. Treatment of the in- 
terstellar medium is crucial for this. Therefore, we look 
to simulations of four star-forming galaxies run with the 
N-Body + smooth particle hydrodynamic code GASO- 
LINE (Wadsley, Stadel, & Quinn 2004; Stinson et al. 
2006) which are described in detail in Christensen et al. 
(2012, 2013) and Munshi et al. (2013). They are able to 
match the sizes, stellar masses, bulge masses, and baryon 
fractions of local star-forming galaxies (Governato et al. 
2007; Brook et al. 2012; Munshi et al. 2013; Christensen 
et al. 2013), and the size evolution of star-forming galax- 
ies since 2 = 1 (Brooks et al. 2011). The simulated galax- 
ies have z = 0 stellar masses which range 4.2 — 4.5 x 10 10 
M 0 for a Kroupa, Tout, & Gilmore (1993) initial mass 
function (IMF). They are referred to as galaxies h239, 
h277, h258, and h285 in the references above. 

These simulations resolve high density peaks compara- 
ble in size to giant molecular clouds. In addition, they 
follow the non-equilibrium formation and destruction of 
H 2 , using a gas-phase and a dust- (and hence metallicity- 
) dependent scheme that traces the Lyman- Werner radia- 
tion field and allows for self-shielding by H 2 gas (Gnedin, 
Tassis, & Kravtsov 2009; Christensen et al. 2013). There- 
fore, the simulations are able to tie star formation to the 
presence of molecular gas, as is indicated by observa- 
tions (e.g., Leroy et al. 2008; Bigiel et al. 2008; Blanc et 
al. 2009; Bigiel et al. 2010; Schruba et al. 2011). The 
efficiency of star formation is tied to the H 2 fraction, 
as described in Christensen et al. (2012). With the in- 
clusion of H 2 -based star formation, stars predominantly 
form at high densities and low temperatures (T < 1000 
K; see also Krumholz, Klein, & McKee 2011; Kuhlen et 
al. 2012). Finally, to take into account the effects of the 
reionization of the universe, an ultraviolet background is 
implemented at 2 = 9, following a modified version of 
the formulation by Haardt & Madau (2001). 

Our simulated galaxies were chosen from lower reso- 
lution simulations to be re-simulated at high resolution 
based on their 2 : = 0 virial masses and because they 
span a representative range of halo spin values and ac- 
cretion histories. The re-simulations follow a 50 Mpc 
co-moving box around the galaxies at lower resolution 
than the galaxy itself. In this manner, cosmological ef- 
fects are fully implemented. The re-simulations are run 
from z = 150 to 2 = 0. 

The spline force softening in the high resolution re- 
gion is 174 parsecs, and is kept fixed in physical parsecs 
since 2 = 10. In the high resolution region, dark mat- 
ter particles have masses of 1.3 x 10 5 M 0 , gas particles 
start at 2.7 x 1O 4 M 0 , and star particles are born with 
30% of the mass of their parent gas particle, which cor- 
responds to a maximum initial mass of 8100 M 0 . Each 
simulated galaxy has ~ 5 million dark matter particles 


within its virial radius at 2 = 0 and more than 14 mil- 
lion total particles (dark matter, gas, and stars). Our 
simulations have the same mass resolution as the Eris 
simulation (Guedes et al. 2011), but include metal line 
cooling (Shen, Wadsley, & Stinson 2010) and the physics 
of molecular hydrogen (Christensen et al. 2012). The 
simulations are run for a WMAP 3 Year Cosmology: 
f! m = 0.24, A = 0.76, H 0 = 73 km s" 1 , ct 8 = 0.77 
(Spergel et al. 2007). 

2.1. Warm & Cold Gas in the Simulations 

The nebular emission lines in real galaxies, from which 
internal kinematics are measured, largely come from H II 
regions, with a small contribution from diffuse ionized 
gas. Gas in H II regions has been ionized by young stars 
after the stars and gas emerge from their dusty birth en- 
vironments inside dense regions of molecular clouds. We 
investigate two types of star forming gas in the simu- 
lations which should bracket the kinematic behavior of 
H II regions: (1) cold gas in the galaxy disks which we 
refer to as “cold gas”, and (2) dense gas in the galaxy 
disks which has been heated by supernovae feedback and 
which we refer to as “warm gas.” Cold gas is defined as 
having temperatures T < 1000K. It has typical densities 
of > 1 amu cm -3 and traces both neutral gas in the disk 
and molecular clouds at p > 100 amu cm -3 (Christensen 
et al. 2013). Warm gas is defined as having p > 0.001 
amu cm' 3 ); it has temperatures in excess of 20,000K due 
to heating by feedback. Star-forming regions in observed 
galaxies likely contain ionized gas from both the cold and 
warm components because they are being observed after 
the young stars are visible and some feedback has had 
a chance to occur. Therefore, we expect the behavior of 
the cold and warm gas in the simulations to bracket that 
of HII regions. 

3. KINEMATIC EVOLUTION OF SIMULATED GALAXIES 

For each of the 4 simulated galaxies, we measure <r g and 
V ro t at discrete redshifts sampling the range of observa- 
tions in Kassin et al. (2012). We measure the intrinsic 
values of these quantities which are unaffected by obser- 
vational effects such as seeing, slit width, and pixel scale. 
The quantity a g is measured as the average integrated 
velocity dispersion of the gas in the direction perpendic- 
ular to the disk, similar to observations. The rotation 
velocity V ro t is measured as the maximum line-of-sight 
velocity of the gas, also similar to the observations which 
use the rotation velocity on the flat part of the rotation 
curve. We use the actual gas particle velocities rather 
than circular velocities since this provides the most di- 
rect comparison to observations. Rotation velocities are 
all corrected for inclination. 

Median values of a g and V rot at discrete redshifts 
for the warm and cold gas in our 4 simulated galax- 
ies are shown in the top panels of Figure 1. Values of 
these quantities for the individual simulated galaxies are 
shown in the bottom panels. In the top panels, the me- 
dian values are compared to those for an observed mass- 
limited sample of 270 star-forming galaxies from Kassin 
et al. (2012). Qualitatively, the simulated galaxies fol- 
low the same trends as the observations: they increase 
in a g and decrease in V ro t with increasing redshift over 
0.1 < 2 < 1.2. In other words, both the simulations and 
the observations decrease in a g and increase in V ro t with 
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Fig. 1. Top: The redshift evolution of the median integrated velocity dispersion <r g and rotation velocity V ro t of warm and cold gas in 
the simulated galaxies (red and blue points, respectively) is compared with that for an observational mass-limited sample of 270 galaxies 
from Kassin et al. 2012 (grey points). The simulations follow the same disk settling trends as the observations: they increase in cr g and 
decrease in V ro t with redshift (i.e., decrease in o g and increase in V ro t with time). The differences in normalizations between the simulations 
and observations can be attributed to differences in galaxy stellar masses for V ro t and temperatures/regions probed for cr g (see text). Solid 
lines show linear fits given in the text. Error bars on the simulations show the rms scatter. As in Kassin et al. 2012, error bars on the 
observations are calculated by bootstrap re-sampling the data in each redshift bin, since their distributions are non-Gaussian. Bottom: The 
evolution of the warm and cold gas (red and blue lines, respectively) in the four individual simulated galaxies is shown: h285 (solid), li239 
(dotted), h258 (short dashed), and h277 (long dashed). Although the simulated galaxies show the same trends as the observations in the 
median, individually they show significant variation with time, consistent with the large scatter of the observations. 


time over the last ~ 8 billion years, demonstrating that 
the simulated galaxies undergo the disk settling found in 
observations. 

As expected, the normalizations of the simulated and 
observed relations differ. The median values of V ro t for 
the simulated galaxies are greater since they are more 
massive on average than the observed galaxies, and V rot 
generally scales with stellar mass for disk galaxies. The 
masses of the simulated galaxies range 4.2 — 4.5 x 10 10 
M 0 versus 6.3 — 50.1 x 10 9 M 0 for the observations. 
(The simulations and observations adopt Kroupa, Tout, 
& Gilmore 1993 and Chabrier 2003 IMFs, respectively, 
which result in consistent stellar masses.) The normal- 
ization of the o g versus z relations for the warm and 
cold gas are higher and lower than the observations, re- 
spectively. In the simulations the warm gas is affected 
by feedback from supernovae, which results in higher o g 
and leads to hotter temperatures compared to observed 
H II regions. By definition, the cold gas is cooler than 
observed H II regions, which results in lower o g . 

3.1. Quantitative Trends of a g andV rot with Redshift 


Linear relations are fit to the median points in Figure 1 
by performing least-squares fits which take into account 
the errors in o g and V ro t (errors in redshift are negligi- 
ble). When fitting, to avoid covariances, we zero-point 
the medians near the middle of the samples such that 
they vary around ~zero. We obtain the following fits for 
cr g for the warm gas, cold gas, and observations, respec- 
tively: 

log a g - 50 = (19.3 ± 6.2 ){z - 0.6) - (4.2 ± 1.8), (1) 

log o g - 20 = (6.4 ± 6.8)(z - 0.6) - (2.3 ± 1.8), (2) 

log o g - 35 = (26.0 ± 5.7 ){z - 0.6) - (3.1 ± 1.4). (3) 

These fits have y 2 values of 1.9, 2.5, and 4.7, respectively. 
As for the observations, the median o g of both the warm 
and cold gas grows with increasing redshift to z = 1.2 
(i.e., decreases with time), although the relation for the 
cold gas is consistent with no evolution. The slopes of 
the relations for the warm gas and observations are con- 
sistent within uncertainties. As mentioned above, the 
normalization of the relations for the warm and cold gas 
are higher and lower than the observations, respectively. 
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z = 0.3 


z = 0.5 


z = 1 



Fig. 2. — Multi-color images of two of our four simulated galaxies are shown at discrete redshifts spanning the redshift range of Figure 1. 
The morphologies become progressively less disturbed with time and the disk grows (from right to left), reflective of the kinematic evolution 
which shows disk settling with time. The images are for galaxies h277 (top panels) and h239 (bottom panels) and are combinations of g , r, 
and z-band images created with the Sunrise radiative transfer code. They are 50 kpc (physical) on a side and the central galaxy in each 
image is viewed at an inclination of 45°. 


We obtain the following fits for V ro t for the warm gas, 
cold gas, and observations, respectively: 

log Vrot— 175 = (— 84.7±11.4)(z— 0.6) — (12.0±5.5), (4) 

log Vrot - 200 = (— 47.6 ±11.5) (z— 0.6) — (0.1 ±3.8), (5) 

log V rot -100 = (—48 ± 11.5)(z — 0.6) — (7.8 ± 3.0). (6) 

These fits have x 2 values of 4.7, 7.4, and 9.2, respec- 
tively. As for the observations, the median V ro t of both 
the warm and cold gas decreases with increasing redshift 
to * = 1.2 (i.e. , increases with time). The simulated 
galaxies have a higher normalization (i.e., have faster 
median rotation velocities) than the real galaxies, as ex- 
pected because they are on average more massive than 
the observed galaxy sample, as discussed above. 

3.2. Kinematics and Morphology 

As is the case for real galaxies over 0.1 < z < 1.2 (e.g., 
Flores et al. 2006; Kassin et al. 2007; Yang et al. 2008; 
Kassin et al. 2012), the kinematics of simulated galax- 
ies are reflected in their morphologies. Images of two of 
our four simulated galaxies are shown in Figure 2 at dis- 
crete redshifts over the redshift range in Figure 1. Both 
galaxies have more disturbed morphologies at z = 1.0, 
and become progressively less disturbed with decreasing 
redshift. At z = 0 they appear as ordered disk galaxies 
with little disturbances or peculiarities. This morpho- 
logical transformation is consistent with the kinematic 
settling they undergo (i.e, decrease in a g and increase in 
V ro t with time). 


The images in Figure 2 are produced by taking into 
account the effects of radiative transfer and dust us- 
ing a Monte-Carlo ray tracing program designed to 
pair with hydrodynamic simulations (Sunrise; (Jonsson 
2006; Jonsson, Groves, & Cox 2010). In short, Sunrise 
creates a spectral energy distribution for each star par- 
ticle in the simulation based on its age and metallicity, 
using the Starburst99 stellar population modeling soft- 
ware (Leitherer et al. 1999). The metallicities of the gas 
particles determine how light from the galaxy is attenu- 
ated by dust, and a constant dust to gas ratio of 0.4 is 
assumed. The resulting images are then convolved with 
the SDSS g , r, and i band filters to produce those in 
Figure 2. 

4. CONCLUSIONS 

For the first time we measure the evolution of the kine- 
matics of simulated star-forming galaxies over the last 
~ 8 billion years, over nearly half of the age of the uni- 
verse. Specifically, we measure the evolution of disor- 
dered motions (a g ) and ordered rotation (V ro t) of gas 
in a suite of four simulated galaxies. Measurements are 
compared with recent observations which show the pro- 
gressive settling of gas in disk galaxies over 0.1 < z < 1.2 
from disordered systems into the ordered disk galaxies 
common today (Kassin et al. 2012). We find that the 
simulated galaxies follow the same trends as the obser- 
vations: they progressively decrease in disordered mo- 
tions ( <r g ) and increase in ordered rotation (V ro t) with 
time. Differences in normalization between the observa- 
tions and simulations can be attributed to differences in 
the average stellar mass of the galaxy samples and gas 
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temperatures/regions probed. 

Reproducing these observations is an important bench- 
mark for any theory or simulation of disk galaxy for- 
mation. The observations trace the internal kinemat- 
ics of star-forming galaxies over a significant period of 
time. Internal kinematics directly dictate how the gas in 
galaxies is arranged, and are strongly affected by physi- 
cal processes such as feedback from star-formation, ma- 
jor/minor mergers, and smooth accretion of baryons. 
Our next step in a future work is to place constraints 


on the processes which have the most direct effect on 
disk settling. 
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